Front dynamics in turbulent media 



A.C. Marti,^ F. Sagues,^ and J.M. Sancho^ 
^ Departament d'Estructura i Constituents de la Materia, 
Universitat de Barcelona, 

Av. Diagonal 647, 
E-08028 Barcelona, Spain. 
^ Departament de Quimica Fisica, 
Universitat de Barcelona, 

Av. Diagonal 647, 
E-08028 Barcelona, Spain. 
(February 5, 2008) 



A study of a stable front propagating in a turbulent 
medium is presented. The front is generated through a 
reaction-diffusion equation, and the turbulent medium is sta- 
tistically modeled using a Langevin equation. Numerical 
simulations indicate the presence of two different dynamical 
regimes. These regimes appear when the turbulent flow either 
wrinkles a still rather sharp propagating interfase or broad- 
ens it. Specific dependences of the propagating velocities on 
stirring intensities appropriate to each case are found and fit- 
ted when possible according to theoretically predicted laws. 
Different turbulent spectra are considered. 



I. INTRODUCTION 

Front propagation has been a problem of great interest 
in a rich variety of nonequilibrium phenomena originat- 
ing pattern structures |^-|^ ■ The most usual and studied 
situation corresponds to a stable (planar) front propagat- 
ing in a quiescent homogeneous and isotropic medium at 
a constant velocity. This situation can be modeled easily 
by a reaction-diffusion equation for the order parame- 
ter. Reaction (nonlinear) terms have to present at least 
two possible steady states: one stable and a second un- 
stable or metastable. Then the stable state propagates 
through the unstable one. The region where the variable 
changes abruptly from one state to the other is called the 
interface or front width, whose size is controlled by the 
diffusion coefficient. These fronts propagate at a velocity 
vq which results from the interplay between the chemi- 
cal time scale (reaction) Tchom and the species diffusivity 
coefficient (diffusion) D: vq ~ (^'''chcm)^^^- These facts 
are well known and have been discussed and applied in 
different fields. 

A more complicated situation is found when consider- 
ing front -like patterns propagating in non-quiescent me- 
dia. In such a case, the additional length and time scales 
introduced by the advecting flow will interact with those 
intrinsic to the reaction-diffusion front dynamics, result- 
ing in new and distinctive propagation regimes. Particu- 



larly interesting in this context is the problem of reaction 
propagation in premixed combustible gases [|. Mainly 
due to its practical relevance in combustion processes, 
the subject has motivated numerous theoretical as well 
as experimental efforts j5|-[l0|]. The generic understand- 
ing coming from these studies is that the turbulent front 
velocity vt is larger than vq, to an extent which depends 
on the intensity and spatiotemporal correlations of the 
turbulent flow. Some other issues, however, such as ve- 
locity quenching effects, role of turbulent spectra remain 
somewhat more open. 

An enlarged perspective has been gained recently when 
experiments on liquid phase reactions [ pl have brought 
out a more simplified scenario as compared to that of 
combustion processes. In such liquid reactions the change 
of density and the increasing of the temperature are neg- 
ligible compared with those occurring in typical com- 
bustion processes p^ ]. In this way experimentalists are 
closer to reproduce the set of simplified assumptions in- 
voked by most of the theoretical models. 

An even more controlled situation concerning liquid 
phase reactions can be envisaged when considering front 
propagation under externally imposed stirring condi- 
tions. This is precisely the context we will address in 
this paper. In essence, stirring can be thought of as a 
sort of controlled way to inject energy into a fluid medium 
through random forces of statistical nature 1 13 ij] . Espe- 
cially suited to reproduce such an scenario are the gener- 
ically called |]l5| synthetic turbulence generating models. 
In particular we will employ here an algorithm to gen- 
erate stochastic, turbulent-like, flows satisfying Langevin 
type equations p6| . Specifically, random flows with zero 
mean velocity and statistically isotropic, homogeneous 
and stationary are created independently of the front 
conditions. 

The closest experimental reference we are aware con- 
cerning this question is the recently reported experimen- 
tal work by Ronney et al. dealing with different 
stirred media. From the theoretical side however there 
is no specific literature focusing on stirring effects, so we 
rely on the general references of front propagation under 
turbulent flows. Most of the existing theoretical mod- 
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els refer to the so-called Huygens mechanism valid for 
thin front conditions ("flamelet regime"). Among them 
we mention those based on the so-called G-equation [7|, 
nonlinear averaging of reaction-diffusion equations 
and specially, for the purpose of comparison with our 
results here, the approaches of Kernstein Q , Yakhot 
and Pochcau These authors based on different as- 
sumptions proposed scaling relations for vt in terms of 
the intensity of the flow Uq- Using the so-called flame 
propagation equation, a variant of the well-known KPZ- 
equation, Kernstein and Ashurst Q obtained a relation 
valid for weak stirring intensities and large correlation 
times of the flow. Applying renormalization procedures 
to the G-equation, Yakhot |jl^ found a relation without 
adjustable parameters. Pocheau claimed that both Kern- 
stein and Yakhot laws converge into a quadratic relation 
when invariance over all the scales is assumed . 

On the other hand, when relying on computer sim- 
ulations [p^ , numerical accuracy and finite size effects 
commonly restricts the examined scenarios to small tur- 
bulent intensities. It is also worth mentioning here the 
alternative stochastic method reviewed by Pope In 
such a formulation, a probability density functional is 
introduced and, after invoking some closure approxima- 
tion and elimination of variables procedure, a tractable 
set of partial differential equations is obtained. The most 
important difficulties of this method remain at the fun- 
damental level of the closure approximation. 

In comparison with these generic procedures, what 
we gain with our random flow generating algorithm is 
an easy control over the most representative statistical 
quantities of the advecting flow: its energy spectrum 
and spatiotemporal correlations. This in turn enables 
us to explore the role of these parameters in determin- 
ing the different propagation regimes and particularly 
in the enhancement of the front speed under turbulent 
conditions. In particular, the two basic experiment ally 
identified modes of front propagation under stirring [[12| , 
i.e. the previously mentioned Huygens propagation mode 
and the distributed reaction zone regime, are qualita- 
tively identified and quantitatively analyzed in our simu- 
lations. Moreover we are able to explore for each regime 
the role of a pair of different energy spectra displaying a 
quite different range of energy containing modes. 

The paper is organized as follows. In the next sec- 
tion we summarize the method to simulate the turbulent 
medium. In Sec. [II we introduce the theoretical model 
of front propagation. Numerical results are presented 
and commented according the existing theoretical pre- 
dictions. We devote Sec. to draw some conclusions 
and perspectives. 



II. TURBULENT MEDIA MODELIZATION 

Here we will present a brief summary of the theoreti- 
cal scheme to simulate a statistically homogeneous, sta- 
tionary and isotropic turbulent flow characterized by its 
energy spectrum. As we want to study fronts dynamics 
under different turbulent conditions we will generate two 
types of spectra. 

As a first example we adopt the Kraichnan's (K) spec- 
trum describing a distributed band of excitations 
around a well pronounced peak centered at some well- 
defined wave number ko, 
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In order to take into account the influence of a broader 
spectrum of modes we will modelize another type of tur- 
bulent medium with the choice of the Karman-Obukhov's 



(KO) spectrum which was introduced |20| to study Kol- 
gomorov turbulence. To this end we select from the fam- 
ily of KO spectra, the following form: 
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which has its maximum also at feg. 

As we will assume that the characteristics of the flow 
are not going to be affected by front dynamics, then our 
method will simply reproduce the above mentioned spec- 
tra in a statistical way (see Ref. for more theoretical 
details) . 

Our starting point is a generalized Ornstein-Uhlenbeck 
like Langcvin equation for the stream function ri{r,t), 
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dt 
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where v is the kinematic viscosity and C(r,t) is a Gaus- 
sian white noise process of zero mean and correlation 

{C{viM)C{^2M))^'ieov5{h~t2)5{v^-V2W^. (4) 

According to these definitions, eq and A arc control pa- 
rameters respectively related to the intensity and charac- 
teristic length of the random fiow. Q[A^V^] plays a rele- 
vant role in our approach because it introduces through 
the random stirring forces the spectra we want to use in 
each particular case. In particular for the K spectrum we 
take 
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where A = (3/2) ^/^A:-!. On the other hand, the KO 
spectrum is reproduced through 
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where A = {%lbY'^k^^. 

From the stream function, the 2-d incompressible ad- 
vecting flow is obtained as usual, 



v(r,t)= - 
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This two dimensional version of our algorithm is cho- 
sen for the simplicity in dealing with the numerical sim- 
ulations, but at the same time expresses our belief that 
the basic trends of front propagation in turbulent media, 
that we want to analyze here, are well reproduced in this 
reduced dimensionality. Moreover, recent experiments in 
chemical reactions in quasi 2-d geometries and theo- 
retical analysis [ p^ support also this assumption. 

The energy spectrum E{k) can be expressed in terms 
of the Fourier transformed operator Q{—X^k^), 



(8) 



The stirring intensity, Uq, (uq = dkE{k)) and in- 
tegral time to and length Iq scales of the random flows 
are easily expressed in terms of the functional and para- 
metric specifications for the noise terms left in (H) and 



In particular for the K spectrum, we have 
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And for KO spectrum, we get for the three basic param- 
eters 
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In Fig. |l| three different spectra with the same intensity 
(kinetic energy) are presented: one of Kraichnan type 
(a) and two of the Karman-Obukhov spectrum b) and 
c). We see how Kraichnan' s energy distribution shows 
a pronounced maximum whereas the other spectra are 
broader and show much longer tails. Spectra (a) and 
(c) have the same fco, but (a) and (b) have the same Iq. 
These facts will be of relevance later on to explain the 
dynamical properties of a front propagating in random 
media represented by each one of those spectra. 

For the sake of comparison, we have considered two 
additional and somewhat related, turbulent conditions. 
The first one, hereafter referred on as frozen stirring, 
corresponds to a frozen configuration of the random flow 
generated according the rules above introduced. The sec- 
ond one referred on what follows as periodic flow, rep- 
resent nothing but a periodic set of fixed eddies. It is 
constructed from the single mode stream function 



Vix, y) = Tjo cos( — ) cos( — ), 



(12) 



representing a periodic array of n x n eddies, where N is 
defined from the system size L a,s, N = L/Ax . 

A detailed presentation of the way the algorithm just 
proposed is implemented to simulate turbulent flows can 
be found in 111™. 



III. FRONT DYNAMICS 
A. The model and theoretical basis 

In this section we present the model used to repro- 
duce a chemical stable front propagating in random me- 
dia with prescribed statistical properties. 

The basic ingredient is a reaction-diffusion scheme 
which exhibits two steady states of different stability and 
connected by a front which otherwise propagates in a sta- 
ble (planar) way. There are many possible choices in the 
literature [||J^ and we will take a very simple one which 
corresponds to that of a generic reaction-diffusion equa- 
tion for the scalar field ■0(r, t), used as a relevant variable: 



(13) 



where the nonlinear function /(ip) has a minimum of two 
zeros corresponding to the steady states. In our numer- 
ical simulations we have chosen the nonlinear function 
/{ip) — ifp' — In this case a stable planar front prop- 
agates the stable state ip ~ 1 ("products") into the in- 
vaded metastable one ip = ( "reactants" ) . Two impor- 
tant parameters are the dimensionless propagation rate 
and front thickness which respectively read: 



The front profile has then the form. 
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Our numerical results have been checked in relation with 
different choices of the reaction term. Confirming our 
belief that within our formulation, this choice should not 
be a relevant point, we have explicitly considered the case 
f{ip) = -0(1 — '0)(V' — 0.25) and the same quantitative 
results have been obtained. 

A turbulent flow is superposed now through a convec- 
tive term added to Eq. (O) 



^ = DV^V + V'^ - V'^ - V.[v(r, t)^}]. 



(16) 



As illustrated in Figs. |^a and ||b, numerical simulations 
of Eq. (16) above, reproduce two distinct limiting regimes 
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of front propagation. In these two figures the flow pattern 
is also plotted. Small arrows indicate the local velocity 
vector of the flow. 

In Fig. Ija we see a front in a situation where the typ- 
ical length scale of the flow is larger than the intrinsic 
one associated to the reaction-diffusion dynamics Sq. It 
is clear that the mean size of the observed eddies is larger 
than the front width. Then a distorted front propagates 
with a larger velocity than in the deterministic case but 
maintaining a still rather sharp and well defined inter- 
face. Such a propagation mechanism is known in the 
combustion literature as the "thin flame," "flamelet" or 
"reaction sheet" regime 

The other limit corresponds to a situation for which Iq 
is smaller than Sq. Now the mean size of the eddies are 
comparable or smaller than the interfacial width. This 
is the situation in Fig. ^ where the front is broader 
and more diffused than in the former case. This situa- 
tion is referred in the literature as a "distributed reaction 
zone" (DRZ) regime Also in this case we also observe 
that the front velocity is larger than in quiescent media. 

We want to emphasize at this point that both phe- 
nomenologies are well known in experiments dealing with 
front propagation for isothermal chemical reactions |12|. 

A specific mechanism of front propagation applies to 
each one of the previously identified modes. The com- 
mon rationale behind the "thin flame" mode is based on 
a HP-like argument: the front has the same local struc- 
ture as in the planar case with normal velocity given by 
vq, but its length increases due to wrinkling and this re- 
sults on faster propagation velocities. This can be under- 
stood as a geometrical consequence of the propagation of 
curved interface with a local velocity vq • Denoting re- 
spectively by Lt and Lq the front length and the lateral 
system size, we have 
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On the other hand if we assume that the effect of the 
flow velocity in the DRZ regime is completely reproduced 
by increasing the diffusive transport inside the broadened 
front then such an effect can be incorporated as a renor- 
malization of the diffusion coefficient. In the DRZ mode 
we simply adapt the first fundamental result of Eq. |lj, 
to obtain 



vt\ { Dt 
voj ^ \ D 



(18) 



where Dt is the effective turbulent diffusion. The 
next and most involved step consists, however, in using 
Eqs. (|l^) and (^8|) above to make detailed predictions for 
vt, or its dimensionless form S, as a function of the stir- 
ring intensity uo, or its dimensionless value uq/vq = Q. 
Let us discuss on what follows our numerical results for 
the HP and DRZ regimes. 



B. Results for HP versus DRZ propagation modes 

In Figs. (|^)-(|^) we present our numerical results cor- 
responding to these two modes of front propagation. If 
it is not explicitly indicated the turbulent flow is that of 
K spectrum. On what refers to the HP mode, our first 
task was to check relation (^^. The collected data for 
the different values of Uq are summarized in Fig. |3[ For 
the sake of comparison, we include in this figure results 
obtained for the other two additional stirring conditions 
introduced in Sec. |ll[ the frozen stirring and the periodic 
flow. According to this figure, the geometric argument 
leading to ( p^ seems well-supported by our simulations 
with the simple exception of those situations involving 
very intense periodic flows. Actually, under these last 
conditions the interface is largely perturbed by the pres- 
ence of overhangs whose dynamics contribute positively 
to the computed velocity, measured as the time varia- 
tion of the rate of occupation of the -0 = 1 state, but 
negatively to the front length. 

In Fig. ^, results for S{Q) are specifically plotted cor- 
responding to three of the previously mentioned stir- 
ring modes. Figure ^a shows the simulation results of 
5^ versus for a turbulent flow for two different val- 
ues of tf). The theoretical predictions of Yakhot [0, 
(5 = exp(gVS'2)) and Pocheau {S^ ^ I + 2Q^) [0 are 
also plotted. Numerical data of the case of frozen stirring 
are included for comparison. Actually, the numerical re- 
sults fit reasonably well the linear relation S'^ = 1 + 
for the whole range of Q. Nevertheless, the slope of such 
a linear law depends onto, a fact not considered in theo- 
retical predictions of Yakhot and Pocheau. Remarkably, 
these theoretical predictions fit better with our numeri- 
cal results for frozen stirring (limit of very large to)- This 
agreement is not surprising because in Pocheau's analy- 
sis, both theoretical and experimental, a very large value 

for this last flow. 



(17) of to was considered |10 



Moreover, 

the behavior largely depends on the examined range of 
Q: the former behavior transforms into a S* = 1 + a'Q's 
when approaching the smallest values of Q here consid- 
ered (Fig. |b). 

Finally, results for the periodic stirring also show a 
crossover from a linear dependence at large Q towards 
the quadratic form: S'^ = 1 -|- a"Q'^ at small Q (Fig. ^). 
Results for the last two cases are in agreement with the- 
oretical predictions by Kernstein et al. 

A final comment is worth emphasizing at this point: 
Turbulent propagation velocities fall always lower than 
those obtained for frozen and periodic stirring. This is 
somewhat at odds with what is reported in Refs []6|,p^, 
although one should keep in mind that the range of Q 
values there considered were of two orders of magnitude 
larger than ours. 

In Fig. I results for the DRZ regimen are presented. 
The numerical data obtained under different stirring 
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mechanisms are plotted. These numerical results are 
compared with the corresponding theoretically predicted 
values based on Eq. (^8|) above. The values of the ef- 
fective diffusion coefficient Dt have been obtained in- 
dependently from direct simulation of pure scalar diffu- 
sion (without reaction). Theoretical predictions based in 
Eq. ( p^ ) exhibit a remarkable agreement with numerical 
results for a broad set of the Q values here considered and 
irrespective of the type of stirring flow considered. More- 
over, the theoretical dependences of Dt on Mq either for 
random flows in the weak stirring limit (Dt — D ~ Mpto) 
p6[ or for periodic flows in the limit of small Peclet num- 
ber {Dt ~ D ^ ^o/D) are clearly confirmed in this 
Fig.|. 

C. Results for K versus KO spectra 

In order to investigate the influence of the differ- 
ent spectra on the enhancement of the front speed, we 
present in Fig. || numerical results corresponding to the 
two regimes of front propagation previously identified 
and the pair of turbulent spectra introduced in Sec. ||. 
Respectively, Fig. |^a reproduces conditions of thin front 
propagation, whereas Fig. 0b those for distributed reac- 
tion fronts. Both in Figs.^a and circles stand for 
results obtained with K spectrum, whereas squares and 
romboids correspond to the KO distribution. Since both 
spectra can be compared either on the basis of their in- 
tegral length scale, Iq, or of their maximum, kg, we have 
chosen to separately consider both cases. Specifically, 
squares represent simulations from both spectra with the 
same value of Iq, i. e. comparing cases a) and b) in Fig. ||, 
and romboids correspond to both spectra taken at equal 
values of ko, i.e. comparing cases a) and c) in Fig. 0. 
On the other hand, black circles in Fig. ^ and black 
and white ones in Fig. are respectively replottcd from 
Figs. ^ and 1]. Just for the sake of comparison, we have 
also included in Fig. ^ as white circles, results for ob- 
tained from the simulated values of the effective diffusion 
through the use of expression (|l^) . 

Let us start with the conditions of the HP regime 
(Fig. ||a). When both spectra are taken at equal Iq, we 
find smaller values of the turbulent front velocity for KO 
than for K spectrum. The situation is just the reverse 
one when both spectra are taken with the same maxi- 
mum. Both findings admit a rather direct interpretation. 
From Fig. ^a we conclude that when both spectra have 
equal integral length scale, the KO distribution displays 
its maximum at a smaller value of fco- This means that in 
this last case larger length scale modulations are going 
to be the most relevant ones in bending the propagat- 
ing interface. This in turn means that the front length 
is going to be less enhanced and, since we are referring 
to the HP propagation mode, so will happen with the 
effective front speed. Contrarily when the two maxima 



coincide, what makes the difference is the long tail in the 
KO distribution. Such large wavenumbers, small spatial 
modulations, are going to be very effective both in wrin- 
kling more intensively the front interface but even more, 
by interfering the front dynamics at scales comparable 
to the front thickness. Both effects lead to larger val- 
ues of the turbulent front velocity. In particular the last 
mentioned feature is evidenced when observing that re- 
sults for the KO spectrum approach in this case those 
that would correspond to a front subjected to the K dis- 
tribution and propagating under the enhanced diffusion 
mechanism (white circles). A somewhat similar enhance- 
ment of the front velocity at moderate Q has been found 
in a theoretical analysis when a thin front is slightly 
perturbed at scales smaller than Sq. 

Let us turn now to the DRZ regime of Fig. ||b. What 
we observe in this situation is that no appreciable effect is 
found when both maxima coincide in the respective spec- 
tra. However, smaller values of S*^ are observed for equal 
integral length scales. Analogously as before, both re- 
sults can be interpreted just by looking at Fig. |[ In the 
first case, one can imagine that the interface thickness 
would be wide enough to comprise the most energetic 
wavenumbers of both flows in such a way that the effec- 
tive diffusion would be similarly renormalized. Contrarily 
when prescribing an equal value of Iq, those more ener- 
getic modes in the KO distribution are shifted to smaller 
wavenumbers and in turn are going to be less effective in 
enhancing the effective diffusion. 

IV. SUMMARY AND PERSPECTIVES 

In summary, we have presented here a simple model of 
a chemical front propagating in a turbulent media. Al- 
though artificial, the model reproduces realistic scenarios 
of front propagation modes in chemical experiments. The 
observed dependences of the propagating front velocities 
on the parameters and characteristics of the flow com- 
pares well with some of the most well-known relations 
proposed in the literature. The simplicity and versatility 
of the computer implementation of turbulent flows and 
the simple model used to describe the front propagation 
in these media open new perspectives in the study of 
more complicated situations in turbulent fluids. 
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FIG. 1. Comparison of Kraichnan 

(a) and Karman-Obukhov (b) and (c) spectra. Continuous 
lines correspond to the analytical expression Eq. (^) and sym- 
bols stand for simulation results in the discrete lattice, (a) 
Kraichnan, Iq — 4.0, ko — 0.27; (b) and (c) Karman-Obukhov 
lo = 4.0, fco = 0.12 and lo = 1.85, ko = 0.27 respectively. 
The intensity (kinetic energy) is fixed to Uq = 0.25. Unless 
otherwise stated a square lattice of 128x128 points and unit 
spacing A = 0.5 has been employed in all our simulations. 
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FIG. 2. Density plot of the stirred reaction front for the K 
spectrum. The arrows stand for the local velocity field plot- 
ted every 4 mesh points in each direction, (a) corresponds to 
the thin front mode, parameter values are D = 0.3, Wo = 3.0, 
A = 4.5 and to = 1.0 (b) corresponds to the distributed reac- 
tion front regime for D = 2.0, ul = 10.0, A = 2.2, to = 1.0. 
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FIG. 3. S vs Lt/Lo for turbulent flow (circles, lo = 4.0, 
to = 0.1), random frozen (squares, lo = 4.0) and periodic ed- 
dies (triangles, 8x8 eddies). In all the simulations presented 
in this figure we have employed D = 0.3. 
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(C) Q 

FIG. 4. S{Q) for three stirring modes: turbulent flows (cir- 
cles, to — 0.1; romboids, to = 3.0), frozen stirring (squares), 
and periodic eddies (triangles). The same parameter values 
of Fig. ^ have been employed. Continuous lines stand for 
the theoretical predictions (see the text), a): Turbulent and 
frozen flows, b): Frozen flow and c): Periodic flow. 
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FIG. 5. S vs Q for DRZ regime, full symbols denote simula- 
tion results of the front propagation while open symbols stand 
for simulation of the pure effective diffusion Eq. ^ (without 
reaction). In all the simulations we have employed D = 2.0. 
Circles, squares and triangles correspond to turbulent flow 
(Zo ~ 2.0, to = 0.1), frozen stirring (Zq = 2.0) and periodic 
eddies respectively (16 x 16 eddies). Continuous lines stand 
for theoretical predictions (see text). 
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FIG. 6. Comparison between K and KO spectra. Circles 
stand for the K spectrum while squares and romboids stand 
for the KO spectrum with the same Zo and the same ko respec- 
tively (see Fig. |^) . Black circles are the same simulation data 
of Figs, ^a and 5, and white circles correspond to the sim- 
ulated effective diffusion, (a) Thin front propagating mode, 
(b) Distributed front propagating mode. For KO spectrum 
we have employed TV = 256. 
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